library(ggpubr)
library(tidyHeatmap)
library(ggrepel)
library(tidyverse)

# abundace barplot ---------
abs_abu2 <- read_tsv('mission/Duan_GSDMx/qiiall/abundance_tables/abs-abund-table-4.tsv',
         skip = 1, name_repair = 'universal')

abs_abu2_tidy <- abs_abu2 |>
  pivot_longer(-1, names_to = 'sample', values_to = 'abundance') |>
  mutate(phylum = str_remove(.OTU.ID, '.+;') |>
           fct_lump_n(n = 10, w = abundance) |> fct_reorder(abundance))

abs_abu2_tidy$phylum |> table()

abs_abu2_tidy |>
  mutate(batch = str_extract(sample, 'batch.'), group = str_extract(sample, 'KO|WT')) |>
  ggplot(aes(sample, abundance, fill = phylum)) +
  geom_col(position = 'fill') +
  scale_fill_brewer(type = 'qual', palette = 'Paired') +
  theme_pubr(legend = 'right', x.text.angle = 45) +
  facet_wrap(~batch+group, ncol = 4, scales = 'free_x') +
  labs(fill = 'Order')

# alpha diversity ---------
alpha.metric <-
list.files('mission/Duan_GSDMx/qiiall/diversity/alpha_diversity/',
           pattern = 'metadata', recursive = T, full.names = T) |>
  map(read_tsv, comment = '#') |>
  reduce(left_join) |>
  pivot_longer(4:7)

alpha.metric |>
  mutate(batch_group = str_c(batch, group)) |>
  ggplot(aes(batch_group, value, color = batch_group)) +
  geom_boxplot() +
  geom_jitter(width = .1, height = 0) +
  facet_wrap(~name, scales = 'free_y') +
  scale_color_brewer(type = 'qual', palette = 'Paired') +
  theme_pubr() +
  labs(title = 'Alpha diversity', y = 'Index value')

# beta diversity ------
beta.name <-
  list.files('mission/Duan_GSDMx/qiiall/diversity/beta_diversity/',
             pattern = 'raw_data', recursive = T, full.names = T)

beta.metric <-
  beta.name |>
  str_extract('(?<=//).+distance') |>
  set_names(beta.name, nm = _) |>
  map(read_tsv) |>
  list_rbind(names_to = 'metric')

beta.metric |>
  filter(SubjectID1 > SubjectID2) |>
  group_by(metric) |>
  mutate(normalized.dist = (Distance - mean(Distance))/sd(Distance)) |>
  ggplot(aes(SubjectID1, SubjectID2, fill = normalized.dist)) +
  geom_tile() +
  facet_wrap(~metric) +
  scale_fill_viridis_c() +
  theme_pubr(x.text.angle = 90, legend = 'right') +
  labs(x = 'sample1', y = 'sample2', title = 'Beta diversity heatmap')

beta.metric |>
  filter(SubjectID1 > SubjectID2, str_detect(metric, 'bray')) |>
  heatmap(SubjectID1, SubjectID2, Distance)

# ancom volcano ------
ancom.batch6 <-
  read_tsv('mission/Duan_GSDMx/qi2var/ancom/Category-batch-level-6/data.tsv')

ancom.batch6t <-
read_tsv('mission/Duan_GSDMx/qi2var/ancom/Category-batch-level-6/ancom.tsv') |>
  rename(id = ...1) |>
  left_join(ancom.batch6) 

anc.batch6sig <- ancom.batch6t |>
  filter(`Reject null hypothesis`) |>
  mutate(species = str_remove(id, ';+$') |> str_remove('.+;'))

ancom.batch6t |>
  mutate(type = case_when(`Reject null hypothesis` & clr > 0 ~ 'Up in batch1',
                          `Reject null hypothesis` & clr < 0 ~ 'Up in batch2',
                          .default = 'NS')) |>
  ggplot(aes(clr, W, color = type)) +
  geom_point() +
  geom_text_repel(data = anc.batch6sig, aes(label = species), color = 'black') +
  theme_pubr(legend = 'right') +
  scale_color_manual(values = c('grey','blue','red')) +
  labs(title = 'ANCOM differential species taxa: batch2 vs batch1',
       x = 'Centered Log-Ratio', y = 'W statistics',
       color = 'Sigificance')

## by group
ancom.ko6 <-
  read_tsv('mission/Duan_GSDMx/qi2var/ancom/Category-knockout-level-6/data.tsv')

ancom.ko6t <-
  read_tsv('mission/Duan_GSDMx/qi2var/ancom/Category-knockout-level-6/ancom.tsv') |>
  rename(id = ...1) |>
  left_join(ancom.ko6) 

anc.ko6sig <- ancom.ko6t |>
  filter(`Reject null hypothesis`) |>
  mutate(species = str_remove(id, ';+$') |> str_remove('.+;'))

ancom.ko6t |>
  mutate(type = case_when(`Reject null hypothesis` & clr > 0 ~ 'Up in WT',
                          `Reject null hypothesis` & clr < 0 ~ 'Up in KO',
                          .default = 'NS')) |>
  ggplot(aes(-clr, W, color = type)) +
  geom_point() +
  geom_text_repel(data = anc.ko6sig, aes(label = species), color = 'black') +
  theme_pubr(legend = 'right') +
  scale_color_manual(values = c('grey','blue','red')) +
  labs(title = 'ANCOM differential species taxa: WT vs KO',
       x = 'Centered Log-Ratio', y = 'W statistics',
       color = 'Sigificance')
